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We consider the problem of quantum multi-parameter estimation with experimental constraints 
and formulate the solution in terms of a convex optimization. Specifically, we outline an efficient 
method to identify the optimal strategy for estimating multiple unknown parameters of a quantum 
process and apply this method to a realistic example. The example is two electron spin qubits 
coupled through the dipole and exchange interactions with unknown coupling parameters - explicitly, 
the position vector relating the two qubits and the magnitude of the exchange interaction are 
unknown. This coupling Hamiltonian generates a unitary evolution which, when combined with 
arbitrary single-qubit operations, produces a universal set of quantum gates. However, the unknown 
parameters must be known precisely to generate high-fidelity gates. We use the Cramer-Rao bound 
on the variance of a point estimator to construct the optimal series of experiments to estimate 
these free parameters, and present a complete analysis of the optimal experimental configuration. 
Our method of transforming the constrained optimal parameter estimation problem into a convex 
optimization is powerful and widely applicable to other systems. 



I. INTRODUCTION 

The tremendous allure of quantum information pro- 
cessing has fueled recent progress in the experimental 
and theoretical understanding of physical systems oper- 
ating in regimes where classical physics fails to hold. The 
precise control and characterization of physical systems 
demanded by quantum information processors, e.g. for 
performing high-fidelity quantum gates, has extended our 
mastery of optical, gas phase, and condensed phase phys- 
ical systems. 

Typically, as the precision to which one must char- 
acterize a physical system increases, the sophistication 
of the techniques used to study the system must also 
increase. Recent tour-de-force experiments have fully 
characterized quantum systems of small dimension by 
performing exhaustive process tomography (e.g. [TJH]). 
Such exhaustive tomography requires resources that scale 
exponentially with the dimension of the system being 
studied and so is infcasiblc for systems much larger 
than those already characterized in this manner. Con- 
sequently, many techniques for approximate characteri- 
zation of large dimensional quantum systems have been 
formulated in recent years |3J HI [5J [B] . 

In many situations one is not completely ignorant 
about the dynamical system being studied. An exper- 
imentalist may have partial knowledge of the system 
through information from system preparation or prior 
characterization studies. In such cases the system char- 
acterization often becomes a problem of parameter esti- 
mation, and an important question arises: how does one 



design an experiment to identify the unknown param- 
eters of the dynamical process most efficiently, or even 
optimally with respect to some metric? Experiment de- 
sign for optimal parameter estimation in quantum sys- 
tems is a natural extension of the equivalent classical 
design problem; one typically attempts to rapidly reduce 
the variance in the unknown parameters by performing 
as few experiments as possible. The goal of experiment 
design is to identify the input states to probe the dy- 
namical process with and the measurements to perform 
on the outputs, so that the variance in the unknown pa- 
rameters can be decreased as quickly as possible (with 
the number of experiments performed). Analytical and 
numerical methods for optimal experiment design have 
been widely explored for one parameter quantum pro- 
cesses (e.g. E U HQl HU Q2]), but very few analytic 
optimality results exist for the multi-parameter case (for 
exceptions, see Refs. [T31 E]), in part because of diffi- 
culties in optimizing over noncompatibile (noncommut- 
ing) quantum observables |15j . Numerical approaches 
to optimal experiment design for quantum tomography 
(when all parameters of the quantum process or state are 
unknown) and Hamiltonian parameter estimation using 
convex optimization were first proposed in Ref . |16j , and 
applied to experiments in Refs. [T7J [15] . The method 
follows from the optimal experiment design approach de- 
scribed in Ch. 7.5 of Ref. [19]. Recently, a similar nu- 
merical approach to multi-parameter quantum process 
estimation, using convex optimization, was formulated 
[SU] and we shall further refer to it below. Experimen- 
tally motivated techniques for multi-parameter estima- 
tion have also been proposed [22] , but the optimality 
and asymptotic performance of these are unknown. 
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In this paper we examine this problem of optimal 
multi-parameter estimation for quantum processes when 
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there are constraints on the possible input probe states 
and on the possible measurements. The constraints 
on the input states and spin measurements result from 
experimental limitations on the types of input states 
(measurements) that can be realistically prepared (per- 
formed). We consider a concrete example motivated by 
an experimental platform for quantum information pro- 
cessing: donors in semiconductors with electrical control 
and measurement [23J ■ We solve the problem of 

precisely identifying the coupling between two electron 
spin qubits that interact through a combination of ex- 
change and dipole-dipole interactions by a preparation 
of input states and measurement of electron spins after a 
suitable interaction period. Note that precise knowledge 
of the qubit-qubit interaction is crucial for the execution 
of two-qubit gates which typically work by transforming 
this interaction into the desired gate by single qubit ma- 
nipulation pulses |26j . We apply a recently re-formulated 
numerical approach to optimal experiment design for 
multi-parameter quantum estimation |16j which also in- 
corporates available experimental configurations into a 
convex optimization |20j . This formulation allows us to 
efficiently identify the optimal characterization experi- 
ment and estimate the number of experimental runs nec- 
essary to achieve a desired accuracy in the estimated pa- 
rameters. 

In section [IT] we provide background on the quantum 
parameter estimation problem and recap the formulation 
of the multi-parameter constrained estimation problem 
as a convex optimization from Refs. [THlfSUj . Section III 
presents the experiment design framework in full gener- 
ality and sketches an algorithm for optimally estimating 
a set of unknown parameters of a quantum process. Sec- 



tion IV introduces the example we explicitly solve in the 
paper: two coupled electron spin qubits. We summarize 
the experimental capabilities of this implementation of 
quantum information processing and give a detailed de- 
scription of the coupling dynamics. Then in section[V|we 
apply our experiment design framework to formulate the 
optimal estimation scheme for identifying the unknown 
parameters under the given constraints. 



The generalization to biased estimators is well known, 
but needlessly complicates our discussion. 

However, some probability distributions are more eas- 
ily estimable than others. Take for example a Dirac-delta 
distribution centered at xq, so that the probability den- 
sity function is p Xo (x) = 5(x — xq). The parameter to be 
estimated in this case is xq and only a single measurement 
is required. On the other hand, accurately estimating the 
mean of a large-variance Gaussian distribution requires 
many samples. Estimability of a parameter is thus a 
property of the probability distribution and is indepen- 
dent of the estimator used. This idea is encapsulated by 
the Cramer-Rao bound, which places a lower limit on the 
variance of any single-parameter estimator |27j . 



varT e (X) > 



1 
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Here, iV is the number of samples and F(0) is the Fisher 
Information, defined as a functional of the probability 
distribution, 



f(0) ^ n-m M x) 



where we use the shorthand pe(x) for the conditional 
distribution p(x\9) and (f(x)) — J f (x) p{x\6) dx (or 
Si f( x i)p( x i\@) if the probability distribution is discrete) 
is the expectation value of f(x). Note that the Fisher 
information is a function of the true value (not the es- 
timate) of the parameter. Intuitively it represents the 
amount of "information" about the parameter in the con- 
ditional probability distribution for the data. 

In the multi-parameter case the generalized Cramer- 
Rao inequality bounds the covariance matrix of the (now 
vector- valued) estimator [27] . 



cov e T e (X) > 



AQ)- 1 
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where 2(9 ) is the Fisher information matrix: 



II. 



PARAMETER ESTIMATION 



1(9) = ((Velnpe(x)) (V e lnp e (x)) T 



Suppose a sequence of data that is independent and 
identically distributed (iid) is drawn from a distribution 
that is parametrized by one or several unknown quan- 
tities. For instance, the distribution could be Gaussian 
with unknown mean and variance. The parameter es- 
timation problem is to estimate the value (s) of the un- 
known quantities from the sample data. 

A central task of parameter estimation is the construc- 
tion of an estimator, Tg(X), which maps the sampled 
data, X, to an estimate, 9, of the parameters. In what 
follows, we will assume the use of unbiased estimators, 

(0) = (T 9 (X)) = 9. 



We have used the notation Ve/(0) = 

(del/' 5el/' ■ ■ • ' IbZ^) ■ 

The Cramer- Rao inequality provides a bound on how 
well we can do when estimating the parameter(s) from 
the data. While the actual variance in the parameter 
estimate is dependent on the particular estimator used, 
there exist estimators that are known to saturate this 
bound asymptotically (in the limit of large N) [57] . An 
example, that we shall employ below, is the the maximum 
likelihood estimator (MLE). The MLE is defined as 
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(X) = argmaxpe(A") 
9 



(2) 



3 



III. OPTIMAL EXPERIMENT DESIGN FOR 
QUANTUM PARAMETER ESTIMATION 

Up to this point we have discussed the mathematics of 
parameter estimation. The physics of a particular prob- 
lem becomes important only in calculating the probabil- 
ity distributions (and their derivatives). Quantum me- 
chanics provides the tools with which these distributions 
can be obtained. 

We begin by defining an experiment, £, as a choice 
of the initial state, po! evolution time, t; and a positive 
operator valued measure (POVM), M = {Mi} |S8|. The 
POVM, also known as a generalized measurement, sat- 
isfies J^-Mj = 1 and Mj > 0, < i < n out . Each 
Mi corresponds to a possible outcome from applying the 
measurement M. Through the application of the Born 
rule, each experiment determines a parametrized family 
of discrete probability distributions, 



pf(i) = Tr (Mi (Ug(t)p Ug(tf)) 
= Tr (M lPe (t)) 



(3) 



Here U g (t) = T exp (-i J* H g (t')dt' /hj is the unitary 

time evolution operator and Hg(t) is the Hamiltonian 
whose parameters, 8, we wish to estimate. pg(i) is the 
probability, given a fixed experiment £ — {po,M,i} and 
assuming the parameter takes the value 8, that we get 
the measurement result i. From this probability distri- 
bution one can calculate the Fisher information matrix 
associated to this experiment: 



l£[e) = y (VflPf W) (VsgfMZ 



(4) 



Inserting this quantity into Eq. Q gives a lower bound 
on the variance of our estimate. We will restrict our 
discussion to closed-system (i.e. Hamiltonian) evolution 
and, for the sake of clarity, to finite-dimensional Hilbert 
spaces. The generalization to non-unitary processes is 
straightforward; the most difficult step being the calcula- 
tion of Vep, which must often be performed numerically. 

It is often the case that an experimentalist has access to 
a number of different initial conditions and measurement 
bases. We would like to answer the question: which of 
these initial conditions and measurements should the ex- 
perimentalist use in order to best estimate the unknown 
parameters in the quantum process? In order words, we 
would like to design our experiment so that we sample 
the quantum process in a manner that produces the most 
information about the unknown parameters. Formally, 
suppose we are given a menu of possible experiments, 
{£} and each time we sample our quantum process, an 
experiment, £ = {p$ , NL £ , t £ }, is chosen with probability 
Af (so J^e Af = !)• The result of that measurement is 
governed by the probability distribution, pf (•), and so is 
associated with its own Fisher matrix - i.e. Eq. Q. The 
probability of any particular measurement result must 



now be scaled by the probability that a particular exper- 
iment will be performed. So the Fisher matrix for the 
combined experimental scheme defined by £ and Ag is 



m = E 



(V fl A£pf (*)) (V fl A £P g(i)) T 



£,i 



= 5> £ X £ (0) 



It is now natural to ask, given a menu of experi- 
ments, what choice of A minimizes Tr = 

Tr (J2s Xs1 e {8)) 1 and thereby provides the best upper 
bound on the average of the estimate variance across all 
parameters through Eq. Q. This optimization problem, 
known as A-optimal experiment design, can be written as 



minimize Tr 



subject to A, > 0, YjAj = l 

i 

Note that the optimization parameter is the vector of 
probabilities Xg. This optimization is difficult because 
the cost function is not linear or convex in the optimiza- 
tion parameter. However, through the use of the Schur 
complement (see Appendix [A]), it can be reformulated as 
the convex optimization problem: 

minimize Tr Q 

subject to ^ ^)^°' F = E X £ l£ ( e )> 

A 4 >0, 53^ = 1. (5) 



To make this problem tractable, the menu of experi- 
ments can be chosen as a discretization of the continuous 
space of all possible experiments. The exact nature of 
the discretization must be determined for each problem 
individually, but, in general, a finer grained discretization 
produces a larger optimization problem. A coarser grain- 
ing will result in a smaller optimization problem, but one 
whose solution will more poorly approximate the true 
achievable lower bound on the variance given by Eq. ([1]) . 
In practice, one will discretize the space of initial states, 
the space of POVMs, and time. Given n p initial states, 
ft M POVM's, and n t times, we have h — n p x hm x n t ex- 
periments, and thus an optimization vector, Xg, of length 
h. 

This procedure for framing the optimal estimation 
problem as a convex optimization over a discrete space 
space of experiments is extremely powerful. Experimen- 
tal constraints can be used to limit the menu of possible 
experiments and the optimal distribution can be found 
quickly, even for large problems. Such a restriction to 
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exclude unfeasible experiments is very difficult to incor- 
porate into a continuous optimization technique. From 
the convex structure of the optimization, we also gain 
insight into the expected results. By the complementary 
slackness theorem [T!5], we expect only a small subset of 
the possible experiments to contribute to the optimal dis- 
tribution. This expectation is borne out in the example 
presented in section [Vj 

Given this formulation for identifying the optimal ex- 
periment, we now detail the entire optimal experiment 
design procedure: 

1. Guess parameters. We always need an initial es- 

timate of the unknown parameters with which to 
begin. This assumed value of the parameters, 6 p , 
can be based on prior knowledge about the quan- 
tum process, other studies, or even an educated 
guess. 

2. Enumerate possible experiments. The menu of 

possible experiments £ is dictated by experimental 
constraints. 

3. Calculate Fisher matrices. For each experiment 

on the menu from step 2, the probability distribu- 
tion for the outcome data, pf, (i), and associated 
Fisher matrices X £ (0 p ), must be calculated using 
the assumed value of the parameters. 

4. Perform optimization. The optimization speci- 

fied by Eq. ^ must be performed to obtain an op- 
timal probability distribution of experiments, Xg. 

5. Perform experiments. The unknown quantum 

process should be probed with experiments dis- 
tributed according to Xg. That is, if a total of N 
samples are taken, \Xg /N~\ of them should be using 
experiment £. 

6. Estimate parameter(s). Use the collected data to 

estimate the parameters using an estimator of 
choice. This results in the refined parameter es- 
timate, 6 e . If the maximum likelihood estimator is 
used, we can readily form the likelihood function 
since the N samples are independent - the likeli- 
hood function will be a multinomial distribution: 



W (x) a jv!n^nNw) 



(6) 



where ng is the number of times experiment £ was 
performed (ng — [Ag/AT], and nf is the number 
of times result i (corresponding to POVM element 
Mf) is obtained. The 8 that maximizes this likeli- 
hood function is 6 e , the maximum likelihood esti- 
mate. 

7. Repeat if necessary. This procedure can be re- 
peated, with 6 p in step 1 replaced by e from step 
6. The decision of whether or not to repeat the 



procedure can be based on a number of factors: (i) 
experimental resources, (ii) desired accuracy: if 9 e 
is very different from 9 p then repeating the steps 
is likely to be helpful. Such an adaptive procedure 
will converge on the true value of the parameter(s), 
9 tl through repetition. 

We now illustrate this procedure by treating a specific 
example of constrained multi-parameter estimation that 
is very relevant to quantum computing: the identification 
of coupling parameters in a multi-qubit system. 



IV. DIPOLE- AND EXCHANGE-COUPLED 
QUBITS 

Donors in silicon have been of increasing interest 
in the quantum computing community since the sem- 
inal paper by Kane in Ref. [23]. Most donor based 
quantum computing schemes use the spins of electrons 
bound to donors to encode qubits. Single qubit read- 
out for this implementation is an active area of research, 
but electrically detected magnetic resonance techniques 
[231 12U1 13TJI I3TI 13"2"] are showing potential for delivering 
high-quality single qubit measurements. In order to ex- 
ecute high-fidelity quantum gates, accurate knowledge 
of the coupling Hamiltonian between two donor-bound 
electron spins is required. Given exact knowledge of 
the location of the donors in the substrate, this coupling 
could in principle be computed theoretically. However, 
donors in silicon devices are subject to uncertainty in 
location that is only magnified by subsequent annealing 
processes. Hence it is highly likely that it will be neces- 
sary to characterize the qubit couplings for each device 
separately and therefore an efficient (and preferably op- 
timal) method of doing this characterization is highly 
desirable. As we will demonstrate in the next section our 
constrained parameter estimation scheme is well suited to 
this task because it is numerically efficient and can han- 
dle realistic experimental constraints. Before applying 
our technique we present some details about the physical 
system. 

Two electrons bound to donors implanted in silicon 
will interact through a combination of the dipole and 
exchange interactions. The spin Hamiltonian governing 
dipole coupling between two qubits is 




+ r Ma(r))(T w 8 (2) 



where 7, is proportionality factor relating the magnetic 
dipole moment operator to the Pauli matrices: 

r = r-2 — r% is the vector connecting the two qubits, and 
(O) = O \9) is the expectation value of O over the 
two electron spatial wavefunction, ^(ri,^). 
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The exchange Hamiltonian, a consequence of the 
Coulomb interaction applied to identical spin-^ particles, 
is 

Here, J is the magnitude of the exchange interaction, 
calculable from the localized, single-qubit wavefunctions, 
(j>,ip, by: 

J = e2 /•/■ ^(r 1 )^(r 2 )0(r 2 )V,(r 1 ) dridra 



explicitly, 



The qubits are also subjected to a magnetic field, B 
B z, leading to the Hamiltonian, 



-71 



B a^ - l2 B af\ 



The dipole moments, 7$, of the two qubits may be differ- 
ent due to local inhomogeneities in the substrate or the 
magnetic field. This results in each qubit, even with- 
out any dipole or exchange interaction, having a dis- 
tinct resonance frequency u>i — jiBo/h, with difference 
Ao; = UJ2 — ^1 « uJi,oj2- 

The static interaction, H 0l is presumed to be much 
larger than either Ha or H e , so it is helpful to work in 
the interaction picture (also known as the rotating frame) 
|33j . The effective Hamiltonian in the interaction picture 
is: 



Where, = e^-^V^e* 1 '-^', 



is the 



„-th 



(7) 



Pauli 



= e -Jfa.i>) + e 2ia, a t a (a) 



& («) = _ie- 2 ^*4 a) 



Substituting these into Q, we find many terms propor- 
tional to e ±2l ( UJ i+ UJ 2)t _ These terms are uerg/ rapidly os- 
cillating and will average to zero in a short time. We take 
the rotating wave approximation and neglect these terms, 
keeping only those that rotate no faster than e ± 2lA ^*. 
This leaves us with: 



.(2) 



h.c 



where, 

HF = 2J — 7l72 
TiG = J + 7172 



In the basis {|tt) , III) , lit) > III)}, this Hamiltonian can 
be expressed in matrix form as: 




Hi = h 



( G 

-G 

\ 



o\ 

Fe 2iAot o 

-G 

G J 



matrix in the rotating frame of the or qubit. These are, Hamiltonian is found to be: 



The unitary evolution operator, 

U I (t)=Texp(-if*H(t')dt / /hj generated by this 



Ui(t) = 



-iGt 









_iire-i(^-G)* s in(nt)/n 






-iFe-^-^siniCl^/Ct 
-i(Au-G)t ( cos (ot) + iAuj sm{nt)/n) 








,-iGt 



(8) 



where we have defined ft = \/F 2 + Aj 2 . From this ex- 
pression for the time evolution operator, the Fisher in- 
formation matrices can be computed through Eq. Efy. 



V. OPTIMAL ESTIMATION FOR DIPOLE- 
AND EXCHANGE-COUPLED QUBITS 

In the model described above, the parameters to be es- 
timated are F, G, and Aw. To simplify the presentation 
and for ease of visualization, we assume here that Au 
has been found through standard resonance techniques, 
and focus our attention on the two remaining parameters. 
The general technique is of course valid for any number 
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FIG. 1: (Color online) The single qubit states and measure- 
ments used to probe the process represented on the Bloch 
sphere. 



of parameters (within computational constraints). We 
choose Auj = 1 and take as the true parameter values, 
9 t — (F t = 1.1, Gt = 0.9). We work in units where these 
parameters are dimensionless. 

Realistic experimental constraints for the optimization 
are that the initial states be easy to prepare and the mea- 
surements be experimentally accessible. This is satisfied 
by assuming that all initial states and POVMs are sep- 
arable. Both the initial state and the POVM set can be 
specified by the choice of a Bloch-vector for each of the 
two qubits. To discretize the space of initial states and 
POVMs, the Bloch vectors are chosen from among the 
26 unit-norm vectors, v, of the form 

v = (ax + j3y + yz)j\/ a 2 + (3 2 + 7 2 

where a,(3,-f G {±1,0} are not all zero. These vectors 
are illustrated in Fig. ([I]) 

Given a choice of two Bloch vectors (cr^ ^ and ^Oq 2 ^, 

the density matrix which describes the resulting initial 
state is 



1 



Pa 



r(2] 



Here <r = (a x , cr y , a z ) is a vector formed from the three 
non-trivial Pauli matrices. Similarly, given a choice of 

two Bloch vectors (cr^ \ and ) the corresponding 

POVM elements, which we choose in this case projective 
quantum measurements, are 

M 1 = I(l + ( CT «). CT W)0(l + (^ ) )- (2) ) 



M 2 = - ( I 



r (2) 



r(a: 



M 3 = I(l-(^ ) )- (1) )®(^ + (^ ) )- (2) ) 
M 4 = i(l-(^)-^)0(l-(^)- (2) ) 



These are projectors onto anti-podal points (along the 

axes defined by ^(Tq 1 '^ and (cr^^) on the Bloch spheres 

of the two qubits. The set of initial states and POVMs 
are explicitly enumerated in Appendix [5] For simplicity, 
we will fix the duration of each experiment in the menu 
to t = 1. Therefore, we have n p — 26 2 , um — 13 2 , n t = 1. 

The Fisher matrices are calculated using an initial 
guess Op = (F p — 1, G p — 1) and the optimal experiment 
is then identified using the convex optimization defined in 

n p n M n t = 114244 



III This optimization over h 



section 

experiments takes < 3 minutes on an average, consumer- 
grade desktop computer. The result of this optimization 
is an experimental configuration with only two elements 
of Ag > (X° £ is the probability distribution describing 
the optimal configuration). This means that the optimal 
process probe need only sample using two experimental 
configurations out of the 114244 possible ones. These op- 
timal experimental configurations are shown in Appendix 
[B] The non-zero elements of \° £ are 0.8 and 0.2, which 
means 4/5 of the process probes should be performed 
with one experimental configuration and the remaining 
1/5 with the other. The Fisher information matrix of the 
optimal configuration is: 



1°(0 P ) = 



1.8853 
-0.18431 



-0.18431 
3.3578 



The results of the experiments were simulated using 
t , the actual value of the parameters and the unitary 
transformation given by Eq. (|8| . We sampled the process 
N = 200 times with initial states and POVMs dictated 
by the optimal distribution, A|, and the resulting data 
was used to estimate the parameters using the maximum 
likelihood estimator, Eq. The likelihood function, 

Eq. Q, is plotted for a large range of the parameters 
F and G in Fig. 
surface yields 0° 



2(a) 



Finding the maximum over this 
JF° = 1.10028, G° = 0.8845). The 
estimate of the parameters is extremely close to the real 
values given by t . In addition, we can bound the vari- 
ance of this estimate using the Cramer-Rao bound: 



Var[F e °]+Var[G°]> 



Tr^p)- 1 ) 
200 



= 0.0042 



As noted earlier, we know that as the number of sam- 
ples, N, increases this bound will be saturated. Thus the 
estimation error is controlled and well-known. In figures 
|2(b)| and 2(c)| we show cross sections across the likeli- 
hood function at the estimated values F e and G e . These 
cross sections show how estimation performance is non- 
uniform for F and G. While the value of F is fairly well 
resolved, the likelihood function is highly periodic in G. 
This periodicity reflects the periodic manner in which G 
enters the evolution unitary, Eq. ([8]). We can break the 
periodicity of the likelihood function by varying the time 
for which the quantum process is probed. Figs. 2(d)| - 2(f) 
show that likelihood function and its cross sections when 
the optimal configuration (for t = 1) is used to probe 
the process for times t — 1,1.1,1.4. Again, a total of 
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(a) Log likelihood function attained 
using optimal configuration for a single 
probe time. 




(b) Cross section of the log likelihood 
function in (a) at the value G = G e . 
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(c) Cross section of the log likelihood 
function in (a) at the value F = F e . 
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(d) Log likelihood function attained 
using optimal configuration for multiple 
probe times. 
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(e) Cross section of the log likelihood 
function in (d) at the value G = G e - 
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(f) Cross section of the log likelihood 
function in (d) at the value F = F e . 
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(g) Log likelihood function attained 
using sub-optimal configuration for a 
single probe time. 
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(h) Cross section of the log likelihood 
function in (g) at the value G = G e . 



(i) Cross section of the log likelihood 
function in (g) at the value F = F e . 



FIG. 2: (Color online) The logarithm of the likelihood function, Eq. for a set of simulated data. Darker areas indicate 
a larger likelihood function. Sub-figures (a)-(c) show the likelihood function attained using the optimal experiment design. 
Sub-figures (d)-(f) show the likelihood function attained using the optimal experiment design when the quantum process is 
probed for different times to break the periodicity of the likelihood in G. Figures (g)-(i) show the likelihood function attained 
using a sub-optimal configuration of experiments to probe the quantum process. In sub-figures (a), (d) and (g) the regions of 
white along the F — axis are where the likelihood function is zero and hence its log diverges. 



8 



TV = 200 samples were taken of the process. Periodic- 
ity in G is largely absent in 2(f)| and furthermore, that 
the likelihood function in 2(d) has a dominant central 



peak around the true values of F and G. This technique 
of probing a quantum process for varied times is essen- 
tial when estimating parameters in unitary processes be- 
cause of the potential for parameters to appear in a pe- 
riodic manner in unitary maps. We note that since the 
probe time, t, is actually a parameter of the process it 
should also be optimized over when identifying the op- 
timal experimental configuration. However, we have not 
here included this step in the optimization in the inter- 
est of keeping the search space of the optimization small 
enough to explore within ss 3 minutes on our simulation 
computer. 

To further evaluate the optimal experiment design we 
compare its performance against a sub-optimal estima- 
tion strategy. The sub-optimal strategy we choose is a 
discrete set of initial preparations and measurements all 
aligned along the principal Bloch sphere axes [x, y, z). 
The 12 possible experimental configurations for this sub- 
optimal strategy are listed explicitly in Appendix|S] This 
is a reasonable naive strategy, and we again collected 
-/V = 200 samples with experiments distributed uniformly 
among the 12 possible configurations. The resulting like- 
lihood function is shown in Fig. 2(g) | and cross sections 
of it in Figs. |2(h) and 2(i) Taking the maximum over 
this likelihood surface yields an estimate of the parame- 
ters: s e ° = (F e so = 1.085, G%° = 0.969). This is clearly a 
poorer estimate of the true parameters. We can also cal- 
culate the Fisher information matrix for this suboptimal 
strategy: 



0.5417 0.1662 
0.1662 0.8562 



This Fisher matrix results in the following bound on the 
combined estimation variance: 



VarfF™] + Var[G~ ] > 



Tr^Cflp)- 1 ) 
200 



= 0.016 



The poorer estimate and the larger bound on the vari- 
ance for the suboptimal configuration are clear indica- 
tions of the superiority of the optimal experiment de- 
sign. Furthermore, the number of experimental config- 
urations required to produce a precise estimate of 9 is 
vastly smaller for the optimal design. In Fig. [3j we plot 
the mean squared error of the maximum likelihood es- 
timate as a function of the number of experiments per- 
formed, N. While the MLE for both the optimal and 
sub-optimal configurations approaches the Fisher infor- 
mation bound (provided by the optimal configuration) 
as N — ► oo, the optimal configuration more rapidly ap- 
proaches this bound. Furthermore, the mean squared 
error of the MLE is lower for the optimally configured ex- 
periments for all iV. To achieve the same mean squared 
error, one must perform roughly twice as many experi- 
ments with the suboptimal configuration as are required 



with the optimal configuration for this particular set of 
guessed and actual parameters. 

10 4 MSE 
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FIG. 3: (Color online) Plot of the mean squared error (MSE) 
of the MLE estimator for the optimal (red squares) and sub- 
optimal (yellow diamonds) configurations. Also shown (solid 
blue line) is the Fisher bound for the mean squared error of 
any estimator as given by the optimal experiment. 

To quantify the cstimability of the the parameters in 
this example, we plot the diagonal elements of the inverse 
of the Fisher information matrix as of function of the pa- 
rameters, F and G in Fig. [4] Note that the optimal (con- 
strained) probe configuration has been determined for 
each (F, G) in the plot since this determines the ultimate 
estimation performanc e lim it. Fig. 4(a) shows element 
[1-^niF, G) and Fig. |4(b)| shows element [l^j^F, G). 



As the values of F and G arc changed these elements of 
the inverse Fisher information matrix remain fairly con- 
stant apart from a few notable excursions. This implies 
that the parameters are nearly equally well estimable 
across all possible values. If the optimal probe configura- 
tion is utilized the single sample variance bound is ~ 0.5 
for F and G. If the parameters happen to lie on one 
of the few peaks of [X _1 ]ii(F, G), or [l^^F, G) then 
they arc slightly more difficult to estimate (i.e., a larger 
number of process probes, N, will be necessary to reduce 
the estimate variance) for any possible experimental con- 
figuration. However, note that neither [T~ 1 ]n(F,G) nor 
\I~ 1 ]22{F, G) diverge for any value of (F,G), and hence 
the parameters are always estimable. 



A. Robustness of Estimation Procedure 

Finally, we turn to the issue of the robustness of 
the optimal experimental configuration identified by our 
method. To evaluate robustness, we calculate the inverse 
Fisher information matrix as a function of the parame- 
ters F and G for a fixed experimental configuration (the 
configuration that is optimal for (F,G) = (1,1))- For 
comparison, we also calculate the inverse Fisher informa- 
tion matrix as a function of the parameters for the fixed 
sub-optimal process probe configuration used above. The 
diagonal elements of these matrices are shown in Fig. [5] 
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(a) (1,1) element oil- 1 (F,G) 




In a real experiment, single qubit operations cannot be 
performed perfectly, and as such will always include a 
small error. The state prepared under such a noisy oper- 
ation will be a mixed state that is proximate to the de- 
sired target state. Such inaccuracies in preparation and 
measurement can be easily incorporated into our proce- 
dure by replacing the probe state (POVM measurement) 
constellations with the corresponding achievable mixed 
states (averaged POVMs). Appendix [C| analyzes the spe- 
cific case of small, random gate error in a single-qubit, 
single-parameter estimation problem. The presence of 
such error is shown to increase the number of experi- 
ments required by an amount proportional to a certain 
measure of the error. 




(b) (2,2) element of T~ X (F,G) 

FIG. 4: (Color online) The diagonal elements of the optimal 
inverse Fisher information matrix over a range of values for 
the unknown parameters, F and G. 



These figures clearly show that the optimal configura- 
tion is much more sensitive to parameter variations than 
the sub-optimal configuration. In fact, the single sam- 
ple variance bound for the optimal experiment is quite 
large at some points. This is a consequence of the small 
number of finely tuned experimental configurations uti- 
lized by the optimal experiment. On the other hand, the 
large number of experimental configurations exploited by 
the suboptimal experiment allows for a moderate perfor- 
mance for almost all (F,G). That the optimal experiment 
is rather sensitive to the accuracy of the initial guess em- 
phasizes the importance of going to the adaptive strategy 
mentioned in section III That is, 



as better estimates of 
the parameters are produced, the process probes should 
be adapted to be the optimal configurations for the cur- 
rent guess for the parameters. We expect that this lack of 
robustness of the optimal experiment will be present for 
the vast majority of parameter estimation problems and 
is not a special feature of the example considered here. 
The cost of finely tuning the process probes to optimally 
estimate the parameters based on an initial guess is that 
these probes become less adept at identifying values of 
parameters too far from the initial guess. 

Another important point governing the success of the 
optimization procedure deals with the experimental abil- 
ity to accurately prepare and measure the qubit states. 




(a) [Zf^]" 1 " Mean:2.03 
Min:0.416, Max:27.2 



l n Opt 



Mean:0.981, 



Min:0.300, Max:20.5 





(c) [X" 1 ] Mean:1.32, 
Min:0.234, Max:2.54 



(d) [X^ 2 1 ] Sub Mean:1.30, 
Min:0.213, Max:3.48 



FIG. 5: (Color online) The diagonal elements of the inverse 
Fisher Information matrix over a range of values for the un- 
known parameters, F and G for fixed process probe configu- 
rations. Subfigures (a) and (b) show the diagonal elements of 
the inverse Fisher matrix when the optimal experimental con- 
figuration for the guess (F,G) = (1,1) is used, and subfigures 
(c) and (d) show these matrix elements when the suboptimal 
experimental configuration identified in the text is used. As 
evident from the large deviations in (a) and (b), the small 
number of experiments used in the optimal configuration re- 
duces the robustness of the procedure to errors in the initial 
guess. 



VI. CONCLUSION 

The precise estimation of quantum processes is a key 
ingredient in the engineering of robust quantum infor- 
mation processing devices. For example, to construct 
two-qubit gates for a quantum computer the interaction 
between qubits must be precisely known. This estima- 
tion task is an increasingly demanding one as the scale 
of the quantum process being estimated increases. Thus 
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it is essential to have experimental techniques that use 
minimal resources, but are also accurate. In this work we 
have demonstrated a method for designing the optimal 
experiments for multi-parameter quantum process esti- 
mation. Particular advantages of the method are that it 
can tackle multi-parameter estimation, it naturally incor- 
porates realistic experimental constraints, and that the 
numerical optimization it involves can be implemented 
efficiently. To demonstrate our approach we have applied 
it to the estimation of parameters dictating the coupling 
of two donor electron qubits in silicon. We found the 
optimal experimental configuration among a very large 
candidate set (> 10 5 experiments) and simulated the 
parameter estimation using this optimal configuration. 
The results show that the our method can drastically 
reduce the number of experiments required to perform 
parameter estimation for quantum processes. We com- 
pared the optimal configuration found by our method 
with a sub-optimal approach and quantified the perfor- 
mance improvement of the optimal configuration. We 
also found that while the optimal experiments designed 
by our procedure - which are based on an initial guess 
of the parameters - perform very well, they are very sen- 
sitive to variations in the actual values of the parame- 
ters, and hence lack robustness. However, the general 
algorithm outlined in |III| takes this into account by spec- 
ifying a recipe for adapting the estimation procedure as 
data about the values of the parameters is obtained, and 
hence is capable of compensating for this lack of robust- 
ness in the results of the optimization. 

A useful extension of this work is to investigate the fea- 
sibility of including a robustness measure directly into the 
cost function of the optimization. It remains to be seen 
if this can be done while maintaining the optimization's 
convexity. Robust estimation procedures have been ad- 
dressed in the context of classical control theory [34] and 
the extension of these results to quantum models would 
increase the practicality and appeal of optimal estimation 
in the quantum setting. 

While we illustrated the method here with the example 
of electron qubits in silicon, the general technique of opti- 
mal experiment design for parameter estimation outlined 



in section III is applicable to a wide array of physical sys- 
tems. An interesting avenue for further research would be 
to apply this method to identify the Hamiltonians gov- 
erning small dipole-coupled spin clusters such as those 
probed in recent experiments with diamond |35[ 136] . 

Finally, although the numerical optimization required 
to find the optimal experimental configuration is con- 
vex, and therefore efficient, in the process of applying 
our technique to the example detailed above we noticed 
that current optimization libraries were unable to han- 
dle a very large (h > 150, 000) search space. Therefore 
a possible extension of this work is to use the inherent 
structure in the parameter estimation problem to form a 
smaller optimization program, or possibly to iteratively 
identify the optimal solution. 



APPENDIX A: SCHUR COMPLEMENT 

Consider a nonsingular block matrix, 

A B 



M = 



C D 



where A is an invcrtiblc submatrix. The Schur comple- 
ment of M with respect to A is defined as, 

Mj A = D- CA~ l B, 

A principal theorem in the study of the Schur comple- 
ment [57] says that 



M > 



both A > and M/A > 0. 



Where M > means that M is positive semi-definite. 
Now consider the minimization problem 

minimize Tri* 1-1 
subject to some constraints 

As stated in the main text, this objective function is not 
convex. We then propose a matrix Q > F^ 1 which is an 
upper-bound on the matrix F . This definition implies 
that 

Q-F- 1 =Q- JF _1 J>0. 

And, because F is positive semi-definite, F^ 1 is as well, 
implying Q > 0. So, by the above theorem, 



Q I 
I F 



> 0. 



So we construct the following optimization problem: 

minimize Tr Q 
subject to same constraints, 



Q I 
I F 



> 0. 



The (convex) matrix positivity constraint enforces the 
(non-convex) constraint Q > F , leaving us with a con- 
vex optimization problem (assuming the remaining con- 
straints are also convex). 



APPENDIX B: PROBE CONSTELLATION AND 
OPTIMAL EXPERIMENTS FROM SECTION [V] 

The 26 single qubit initial states available to probe 
the quantum process in the example presented in section 
|V] are given in table |T] Explicitly, for each set of Bloch 
angles, the initial state is: (</>,$)) — cos(^>/2) |0) z + 
e sin(</>/2) \l) z where \i) z are a z eigenstates. The 13 
single qubit POVMs assumed available in section [V] are 
defined by the first 13 angles in table [T] and each POVM 
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*T-Opt 

- L 0.2 



2.03 -0.034 
-0.034 2.82 



has two projector elements given by: \ip((j), 0)) (■0(</>,#)l Fisher matrix, X opt (0 p ) = J2s ^f2^ pt (#) 

wdi-\i>(4>,e)){i>^,ff)\. [ ^ : n 

The sub-optimal estimation strategy used in section 
|V| used a fixed set of state preparations and measure- 
ments to probe the quantum process. There were 12 
possible experimental configurations and they are explic- 
itly enumerated in table [TTJ The initial states and POVM 
elements are defined explicitly in terms of these Bloch 
sphere angles by the same procedure outline above. 



KT 1 = ( 



0.49 0.0059 
0.0059 0.35 



*T-Opt 

X 0.8 



1.85 -0.22 
-0.22 3.49 



TABLE I: Bloch sphere angles for the 26 initial states in sec- 
tionjv] 4> is the polar angle and x = cos~ 1 (l/v / 3). Antipodal 
points are equivalent when choosing POVM's, leading to 13 
inequivalent, single-qubit measurement bases. 



opt] — 1 






e 








tt/4 


{0,tt/2,tt,3tt/2} 


X 


{7r/4,37r/4, 5tt/4, 7tt/4} 


7r/2 


{0, tt/4, tt/2, 3tt/4, tt, 5tt/4, 3tt/2, 7tt/4} 


TT - X 


{7r/4,37r/4,57r/4,77r/4} 


3tt/4 


{0,7r/2,7r,37r/2} 


TT 






TABLE II: Bloch sphere angles (</>, 9) for the 12 experimental 
configurations used by the sub-optimal estimation strategy in 
section [V] (j> is the polar angle, and Ql and Q2 refer to qubit 
1 and qubit 2. 



Init. state Ql 


Init. state Q2 


POVM Ql 


POVM Q2 


(0,0) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 


O,o) 


(0,0) 


(0,0) 


(ir/2,0) 


(-tt/2,0) 


(0,0) 


(0,0) 


(tt/2,0) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 


(0,0) 


(tt/2, tt/2) 


(tt/2, tt/2) 


(0,0) 


O,o) 


(tt/2, tt/2) 


(tt/2, tt/2) 


(tt/2,0) 


(-tt/2,0) 


(tt/2, tt/2) 


(tt/2, tt/2) 


(tt/2,0) 


(0,0) 


(tt/2, tt/2) 


(tt/2, tt/2) 


(0,0) 


(0,0) 


(tt/2,0) 


(tt/2,0) 


(0,0) 


O,o) 


(tt/2,0) 


(tt/2,0) 


(tt/2,0) 


(-tt/2,0) 


(tt/2,0) 


(tt/2,0) 


(tt/2,0) 


(0,0) 


(tt/2,0) 


(tt/2,0) 



The optimal experimental design from the h = 114244 
possible configurations (defined by all possible combina- 
tions of initial state and POVM parameters from Table]!]) 
consists of only two experiments. These are given in Ta- 



ble III and graphical representations of the initial states 
and POVM axes are given in Fig. [6] 

Given below are the Fisher information matrices for 
the experimental configurations chosen by the optimiza- 
tion procedure. Following these is the inverse of the total 



r T opt] 
L x 0.8 J 



0.54 0.035 
0.035 0.29 



[0.2x2^+0.8x2$] '= ( 



0.53 0.029 
0.029 0.30 



TABLE III: Bloch sphere angles ((f), 8) and relative weights in 
A£ for the two experimental configurations that are optimal 
for the estimation problem of section [V] (f> is the polar angle, 
and Ql and Q2 refer respectively to qubit 1 and qubit 2. 





Init. state Ql 


Init. state Q2 


POVM Ql 


POVM Q2 


0.2 
0.8 


(3tt/4, 3tt/2) 
(tt - x, 7tt/4) 


(x,tt/4) 
(x,tt/4) 


O/4,0) 
O/4,0) 


0/4, tt) 
(x,5tt/4) 



APPENDIX C: 



ROBUSTNESS TO GATE 
ERRORS 



When a pure state is acted upon by a noisy gate, the 
result is a mixed state. This mixed state can be repre- 
sented by a Bloch vector which terminates on the inte- 
rior of the Bloch sphere. Though the details depend on 
the error model for the gate, imperfections in preparation 
and measurement can easily be taken into account by our 
formalism. One simply optimizes over a discretized set 
of imperfectly prepared input states and imperfect mea- 
surements. Note that even though all inputs states and 
measurements treated in the example were pure states 
and projective measurement, respectively, our formal- 
ism is not restricted to optimizing over such states and 
POVMs. Specifically, imperfections in state preparation 
can be taken into account by considering the result of 
these imperfections on the Bloch vector, v, of the state, 
defined through it's relation to the single-qubit density 
matrix. 

If we assume that our target state is pure, then vj is 
of unit-norm and the density matrix is: 

P = \ {I + «/ • °) 

Random error in preparation of the initial state corre- 
sponds to the creation of a mixed state. If this error 
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(a) A? = 0.2 experiment 




(b) A| = 0.8 experiment 

FIG. 6: (Color online) Bloch sphere representations of the 
initial states and POVM axes for the two experiments of the 
optimal configuration. The green (dotted) lines are Bloch vec- 
tors for the initial states of each qubit, and the red (solid) lines 
define the axes whose antipodal points define the projectors 
of the optimal POVM for each qubit. 



is assumed to be such that the final state is instead cre- 
ated with some finite probability density surrounding the 
target state, the effect on the Bloch vector is that it is 
contracted by some factor, v 1 * = (1 — e)vf. The details of 
the error govern the magnitude of e. (Of course, there are 
errors which do not just shrink the Bloch vector, but also 
rotate it. As long as these errors are well characterized, 
then a similar analysis may be performed.) The density 
matrix is then. 



p'=\{I+(l-e)v r a) 



In the case where there exists only a single parameter, 9, 
the Fisher information takes the form: 



Pi(6) 



The probability, Pi(6), as given by the Born rule for a 
POVM element, Mi, is 

Pi (6) =Tr (M iP ') 

= ~Tr(M i (I+(l-e)v r a)) 

= ^ (Tr (M^ + (l-e)Tr(v r aM i )) 

Then the Fisher information becomes, in terms of the 
Bloch vector, 



1 = 



E 



Tr (Mi) + (1 - e)Tr (v } ■ aM t ) 



i 

(1 - e) 2 l 



(^(vraMj)) 
Tr (Mi) + Tr (v f ■ ffM t ) 



Here, Xq, is the Fisher Information achieved without the 
presence of gate error. The estimator error, vaxEg, is 
bounded by the Cramer-Rao inequality, 



vaxEg > 



1 



N 



(1 -e) 2 N 



So to achieve the same bound on the estimator variance 
as is found with perfect gates, one must increase the num- 
ber of measurements from N to N' w iV(l + 2e). If there 
are similar POVM errors as well, then a nearly identical 
calculation shows that N' w N(l + 4e). 
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